Critical phenomena in globally coupled excitable elements 
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Critical phenomena in globally coupled excitable elements are studied by focusing on a saddle- 
node bifurcation at the collective level. Critical exponents that characterize divergent fluctuations 
of interspike intervals near the bifurcation are calculated theoretically. The calculated values appear 
to be in good agreement with those determined by numerical experiments. The relevance of our 
results to jamming transitions is also mentioned. 
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The understanding of cooperative phenomena in 
noncquilibrium systems is one of the most important 
problems in physics. In contrast to those in equilibrium 
systems, their nature essentially depends on the dynam- 
ical properties of the systems. This makes it difficult to 
perform a systematic study of the problem. Thus, it is 
necessary to investigate the typical cooperative phenom- 
ena in noncquilibrium systems. 

Recently, critical behaviors have been observed ex- 
jerimentally fil JM EL S S 0| and numerically 
I 111 . Il2l . Il3l . 1 14 Il5l | in typical examples of cou- 
pled excitable elements such as neural networks and car- 
diac tissues. In general, such critical behaviors are clas- 
sified into several groups on the basis of the exponents 
of divergences. The classification enables the realization 
of a universality class for critical phenomena in coupled 
excitable elements. However, the broad distribution of 
the phenomena makes it difficult to elucidate the mecha- 
nism of the criticality. When we consider the role of the 
mean-field Ising model in theories of equilibrium statis- 
tical mechanics, it becomes apparent that it is necessary 
to develop an elegant method for theoretical analysis of a 
minimal model describing critical phenomena in coupled 
excitable elements. 

Thus, with this aim in mind, we analyze a previously 
proposed simple model for globally coupled excitable el- 
ements in this Letter [lf|. In particular, we focus on 
a divergent behavior with respect to parameter change 
around a saddle-node bifurcation because the excitability 
of this model is related to the bifurcation. It should be 
noted that such transition properties have been studied 
for different excitable systems 0,11, HI 11, 3]- The main 
achievement of this Letter is a theoretical derivation of 
the critical exponents that characterize the singular be- 
havior near a saddle-node bifurcation. 

Model: The excitable nature of a system is character- 
ized by the existence of spikes in a time series. Mathe- 
matically speaking, spikes are described by trajectories 
near a homoclinic orbit in a differential equation. As a 
simple example, let us consider an ordinary differential 
equation dt4> = uj — hsm <f> for a phase variable 4> € [0, 2w], 
in which there exists the homoclinic orbit when to = h. 
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FIG. 1: Phase diagram. (L) on the displayed curve satisfies 
0.097 < (L) < 0.107. Inset: (right) Amplitude of oscillation A 
as a function of h. (left) Typical samples of the time evolution 
of Y for h — 1.0 (solid line), h — 1.02 (dashed line), and h = 
1.1 (dotted line). Here, T = 0.05 and N = 100. The arrow 
represents the direction of the parameter change. (L) > 0.1 
in the states where the symbols are placed. 



Then, when h is slightly larger than w, a small perturba- 
tion for the fixed point fa = sm~ (u>/h) yields one spike. 
On the other hand, when h is slightly less than to, the 
system shows an array of spikes with a long interspike 
interval. The qualitative change in the trajectories is an 
example of saddle-node bifurcation. By using this sim- 
ple dynamics as a model of excitable element, we study 
globally coupled excitable elements {fa}^L 1 under the in- 
fluence of noise [Hal : 



K N 

dtfa = oj — h sin fa sin( 

3=1 



(1) 



where £i(t) is Gaussian white noise that satisfies the re- 
lation (&{t)£ 3 (t')) = 2T5 i . J 5(t — t')- Without loss of gen- 
erality we can assume K = 1, and we restrict our in- 
vestigations to the case u> = 1. The control parameters 
are h and T. All numerical results in this Letter have 
been calculated by employing an explicit discretization 
method with a time step St = 0.05. 

The collective behavior of this system is described by 
the time evolution of a complex amplitude, which is given 
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by 



1 N 



j<l>3 



(2) 



In particular, for the expression Z — X + iY, where X 
and Y are real numbers, the expectation values of the 
angular momentum L = X(d t Y) — (d t X)Y are used to 
distinguish the oscillatory states ((L) ^ 0) from the sta- 
tionary states ((L) = 0) In Fig. [TJ we show an 
approximate phase diagram in the form of a curve that 
satisfies the condition 0.097 < (L) < 0.107 in the param- 
eter space (T, h) . A similar phase diagram was obtained 
in Ref. [1 61 ] by measuring the frequency of the time- 
dependent phase distribution function. Here, the curve 
starting from (T, h) ~ (0, 1) is related to a saddle-node 
bifurcation, while the curve from (T,h) ~ (0.5,0) is re- 
lated to a Hopf bifurcation. It should be noted that a 
complicated bifurcation diagram appears near T = 0.1, 
which originates from the Takens-Bogdanov type bifur- 
cation [l8( | . 

Preliminary: In this study, we focus on systems near 
a saddle-node bifurcation. First, we fix T = 0.05 and 
change h across the bifurcation from below. Here, the 
amplitude of oscillation A = ([m&xtX (t) — mintX(i)]) 
changes discontinuously at the bifurcation. (See inset of 
Fig- U) The discontinuous change in the amplitude is 
in a sharp contrast to a super-critical Hopf bifurcation 
at a collective level, where the amplitude of oscillation 
changes continuously at the bifurcation It should 

be noted that in a manner similar to that of the criti- 
cal phenomenon in equilibrium systems, the continuous 
transition leads to a critical diver gen ce of amplitude fluc- 



tuation [20J]. (See also Refs. [2l|, |22j as reviews.) Thus, 
the discontinuous nature of the transition is not indica- 
tive of the appearance of critical phenomena. 

Nevertheless, based on the fact that a typical time scale 
diverges at a saddle-node bifurcation, we take into ac- 
count the fluctuation of intcrspike intervals. Explicitly, 
by using the phase of collective oscillation, 

= arg(Z), (3) 

we define the interspike interval / as the minimum time 
interval [t, t + 1] over which the time integration of dt is 
equal to 2tt for a time t satisfying 0(t) = —tt/2. As the 
most primitive statistical quantities of /, we measured its 
average and fluctuation intensity defined by 



x {h,N) = N [(!')- {I) 



(4) 
(5) 



In order to determine the divergent behaviors near the 
bifurcation in the thermodynamic limit, we performed 
finite-size scaling analysis by using systems with N = 10, 
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FIG. 2: I„N- 1/3 as a function of (1 - h/h c )N 2/3 . Here, 
h c = 1.0283. The guide line represents a power-law function 
with exponent —1/2 
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FIG. 3: X N~ 3/2 as a function of (1 - h/h c )N 2/3 . Here, h c = 
1.0283. The guide line represents a power-law function with 
exponent —5/2 



100, and 1000. For each system, the values of I*(h,N) 
and x{h> N) were calculated for several values of h. Then, 
we assume the scaling relations 



I,(h,N) ~ N$F Z -N^ 
X (h,N) ~ NiF x \ l ^N 



(6) 
(7) 



where the exponents v, £, and 7 and the critical value h c 
are determined so that the scaling relations are valid. We 
also assume that a distribution function of I is expressed 
as a function of IN^/ 1 ' when h = h c . By applying this 
assumption to x(/i, N) in ([5]), we find a relation 7/1/ = 
2£/^ + 1, which yields 



7 -se- 



ts) 



Moreover, since /* and x are independent of N in the 
regime (h c — h)N 1 / v 3> 1, the asymptotic behaviors can 



be derived as Fi{x) 



~ 4 and F x {x) 



With 



the consideration of these conditions, we determine the 
values h c = 1.0283, v = 3/2, C = 1/2, and 7 = 5/2, 
for which the excellent collapses to universal curves are 
found, as displayed in Figs. [2] and [3] 

Theory: We now present a theory for the results 
£ = 1/2 and 7 = 5/2. (y is then determined from ©.) In 
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the argument below, we assume that e = h c — h is a suf- 
ficiently small positive constant and consider the asymp- 
totic limit N — > oo for the assumed value of e. 

We first notice that for a sufficiently small value of 
T, the excitable elements are almost in synchronization. 
Thus, when setting 4n = 9 + 5<fii, we assume that |<5</>,| <C 
1. From this assumption and the definition of 9 given in 
<j3j> , we can derive the equation 



dtO = uj — hi 



17) 



(9) 



with (r](t)r](t')) = 2(T/N)S(t - f), where we have ig- 
nored the contribution of 0(^2^ =1 (S(f>i) 2 /N) to the time 
evolution of 9. Within this approximation. h c is deter- 
mined as h c — to. Although the equation we analyze has 
become quite simple, the calculation of the critical expo- 
nents is still non-trivial. By using a special technique, we 
can derive the distribution function of /, from which we 
can calculate the values of the exponents 23j . Since the 
calculation requires a complicated procedure, we present 
a method by which the values of the exponents can be 
determined without the distribution function of /. 

The basic idea of our analysis is to consider a distri- 
bution function of the average frequency Cl over a time 
interval At = MI*, where M is a large number indepen- 
dent of e. (Note that At depends on e.) For the explicit 
expression 



1 



At 



A/ 



ci = — dt(d t 9), 



(10) 



we can expect a large deviation property, which is given 
as 



p(ci = ci) 



-MNG(n)/T 



(11) 



where the rate function G(Cl) takes a minimum value zero 
when = 0*. Then, it can be shown that I* in (UJ is 
equal to 2-k/CI*. 

We now estimate the rate function G(Cl). Let [6] be 
a trajectory (6(t))£± , and 9(0) is fixed as an arbitrary 
value. The probability density of trajectory is then ex- 
pressed by 



(12) 



where f(9) = ui — h sin 9, the prime represents the deriva- 
tive with respect to 9, and Z is a normalization factor. 
The last term corresponds to a Jacobian term associated 



with the transformation from a noise sequence (r](t))^ 
to the trajectory [9]. By formally expressing P(Cl = Cl) 
as 



p(ct = ci) = / v[0]v([0})6 



dt{d t 9) 



(13) 

we consider the trajectory whose weight becomes most 
dominant in the limit N — > oo. The trajectory, which is 



denoted by 9q, is a periodic solution with period 27r/fi 
of the variational equation d 2 9(t) = —deU(6)/2, where 
U{9) = -f(9) 2 . The solution 9 s n (t) is obtained from the 
energy conservation equation, which leads to the deriva- 
tion of 



8 t 9 s n = JE(fl) - U(9l), 



(14) 



where the parameter E(Q) is related to the frequency CI 
as 



2tt 



2tt 



CI Jo ^E{Cl) - U{9) 



(15) 



Since 9q contributes to P(Ct = CI) much more than 
other 27r/f2-pcriodic trajectories, it is reasonable to ex- 
pect that P(Cl = Cl) ~ V{[9 s n ]). The substitution of Q] 
into (HU) yields 



de WE(Cl) - U{9) y/=U{ff)) 2 



JMcT) - u{9) 



(16) 

It can be observed that dG(Cl)/dCl\Q m = and G(fi*) = 
when Cl* satisfies the condition E(Cl*) = 0. Therefore, 
the rate function G(Cl) takes a quadratic form 



G(Ct) = B(e)Cl- 4 (Cl-Cli,) 2 
when Cl is close to Cl*. where B(e) is calculated as 



B(e) 



8V2tt 



5/2 



+ 0(e 7 / 2 ). 



(17) 



(18) 



Furthermore, by considering (|15[) with E = 0, we obtain 

Cl* = V2~7e 1/2 . (19) 

Now, we consider the average of I during the time 
interval At, which is denoted by J. It can be easily con- 
firmed that J = 2tt/CI. Then, by the transformation of 
the variable in (jTTJ) and pT|) . we derive 



P( J = J) ~ e 



-MNB(e)(J-2n/n*) 2 /(4ir 2 T) 



(20) 



By substituting (JTSJ) and (TTS]) into lf2"0"j). we find that 
(J) ~ e- 1 ' 2 and ((J - (J)) 2 ) ~ e" 5 / 2 . Since these e 
dependences should be equal to those of /* and \, we ar- 
rive at the theoretical results ( = 1/2 and 7 = 5/2. These 
values coincide perfectly with the numerical values. 

Furthermore, our analysis yields a new formula for 
the phase diffusion constant D, which is expressed by 
D = At([Cl- Cl*) 2S J 12 because Cl = (6 (At) - 9(0))/ At. 
Indeed, from pip , we obtain 



D 



T 



2N G"(Cl*) : 



(21) 



which leads to the power-law behavior D = 
(3T/32V2ttN)€- 1 + O(e ). Here, with the crossover rela- 
tion e ~ (N/T)- 2 / 3 , we conjecture D ~ (T/N)(N/T) 2 ^ 3 
at e = 0, which was reported in Ref. [24| . 
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Concluding remarks: We have studied a simple model 
that exhibits critical behavior near a saddle-node bifur- 
cation. The power-law divergence, x — e -5 ^ 2 , which we 
have predicted for coupled excitable elements will be ob- 
served in experimental systems. Complicated systems 
such as those with a tactical network or integrate-and- 
fire dynamics will be analyzed by extending our theory. 

The analysis of finite-dimensional systems is the next 
theoretical problem. As usual in critical phenomena, 
we wish to determine the upper-critical dimension above 
which the values of the exponents are the same as those in 
the globally coupled model. Then, we intend to develop 
a systematic method to take into account non-Gaussian 
fluctuations. The construction of such a theory is ex- 
tremely interesting. 

Before ending this Letter, let us recall that the am- 
plitude of oscillation exhibits a discontinuous transition 
at the saddle-node bifurcation. Here, it should be noted 
that the co-existence of critical fluctuations with a dis- 
continuous transition is one of the remarkable features of 
jamming transitions [25| . This is not an accidental coin- 
cidence and can be explained in the following manner. 

A standard characterization of the critical nature near 
a jamming transition is based on the nonlinear suscepti- 
bility X4(*)i which quantifies the fluctuations of unlocking 
events during a time interval t |26H . Among the several 



theories for Xi(t) 




one theory states that 



the divergence of Xi(t) originates from the critical fluctu- 
ations of the time when an unlocking event occurs [30| . 
By employing the method in Ref. [301 ] . we can discuss 
the divergent behavior of amplitude fluctuations in the 
present problem. That is, the coexistence of a discon- 
tinuous transition and a critical fluctuation in coupled 
excitable elements can be described in a manner similar 
to that in jamming transitions. 

Moreover, it has been recently shown that dynamical 
behaviors of the fc-core percolation in a random graph ex- 
hibit a saddle-node bifurcation at the percolation point 
3l| | . Since it has been known that the fc-core percola- 



tion is related to a kinetically constrained model and a 



random-field Ising model [32, |33j, |3J, |35j, |36| , our work 



might be useful for theoretical analysis of such systems. 

We hope that our theory of the nontrivial behavior of 
a simple model will stimulate further studies on subjects 
that increase the understanding of the cooperative nature 
of nonequilibrium systems. 
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